Cite this Article
Qiao Yao-Bin, Qi Hong, Zhao Fang-Zhou, Ruan Li-Ming. Accurate reconstruction of the optical parameter distribution in participating medium based on the frequency-domain radiative transfer equation. Chinese Physics B, 2016, 25(12): 120201  
Permissions
Accurate reconstruction of the optical parameter distribution in participating medium based on the frequency-domain radiative transfer equation
Qiao Yao-Bin, Qi Hong†, , Zhao Fang-Zhou, Ruan Li-Ming
School of Energy Science and Engineering, Harbin Institute of Technology, Harbin 150001, China

 

† Corresponding author. E-mail: qihong@hit.edu.cn

Project supported by the National Natural Science Foundation of China (Grant No. 51476043), the Major National Scientific Instruments and Equipment Development Special Foundation of China (Grant No. 51327803), and the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (Grant No. 51121004).

Abstract
Abstract

Reconstructing the distribution of optical parameters in the participating medium based on the frequency-domain radiative transfer equation (FD-RTE) to probe the internal structure of the medium is investigated in the present work. The forward model of FD-RTE is solved via the finite volume method (FVM). The regularization term formatted by the generalized Gaussian Markov random field model is used in the objective function to overcome the ill-posed nature of the inverse problem. The multi-start conjugate gradient (MCG) method is employed to search the minimum of the objective function and increase the efficiency of convergence. A modified adjoint differentiation technique using the collimated radiative intensity is developed to calculate the gradient of the objective function with respect to the optical parameters. All simulation results show that the proposed reconstruction algorithm based on FD-RTE can obtain the accurate distributions of absorption and scattering coefficients. The reconstructed images of the scattering coefficient have less errors than those of the absorption coefficient, which indicates the former are more suitable to probing the inner structure.

PACS: 02.30.Zz;45.10.Db;78.20.Ci;78.47.jb
Keyword:inverse problem;adjoint differentiation;multi-start conjugate gradient method;frequency-domain radiative transfer equation
1. Introduction

The theory of inverse problems has been widely developed in the past few decades due partly to its importance of applications in many important scientific and technological fields[1–5] and the advent of the scene of large computers and reliable numerical methods. The reconstruction of optical parameters in the participating medium is a typical inverse radiative problem encountered in retrieving highly-correlated parameters and has attracted significant attention in recent years. For instance, in many fields of engineering and technology, such as nondestructive testing,[6] optical tomography (OT),[7] infrared remote sensing,[8] information processing,[9] and combustion diagnosis,[10] the required optical parameters are unknown and need to be retrieved with the help of the probed radiative signals on the boundary and corresponding inverse theory. Among these applications, the OT has become a research hotspot in the field of the inverse radiative problem due to its unique advantages and wide potential applications. Theoretically, OT is an effective probe technique to reconstruct the distribution of optical parameters by analyzing near-infrared light intensities measured on the boundary of the participating medium. Given that the distribution of optical parameters is closely related to the internal structure of the medium, the optical image reconstructed via OT can also help to explore the internal structure (including the geometric shape, size or location). Compared with the traditional medicine imaging technologies, OT is non-invasive, portable, and capable of producing real-time images of clinically relevant parameters. In the field of biomedical science, the OT technology is considered to be a safer and cheaper supplement to the traditional imaging technologies and has received considerable attention on early tumor diagnosis.[11]

The reconstruction of optical parameters requires considering the laser–tissue interactions to obtain the information contained by the measured radiative intensity. The basic equation to describe the propagation of photons during the laser–tissue interactions is the radiative transfer equation (RTE). Thus an effective forward model needs to be established to simulate the RTE in the participating medium. Theoretically speaking, the forward model applied to OT can be classified into three different techniques according to the incident laser sources: the continuous laser (steady RTE),[12–14] the pulsed laser (time-domain RTE),[15–19] and the modulated laser (frequency-domain RTE).[20–23] Among these three models, the frequency-domain model, which can avoid the technical limitations intrinsic to the use of the time-domain model and gives some additional phase information compared with the steady model, has been regarded as the most promising technique to solve the inverse problems of OT accurately and efficiently. In particular, frequency domain reconstruction can reduce cross-talk between absorption and scattering parameters.[20] In the past few decades, many research groups have conducted considerable researches of the reconstruction of optical parameters based on the frequency-domain model. In the earlier research work,[24,25] the diffusion approximation (DA) model was usually employed as an approximation of the frequency-domain RTE (FD-RTE) to simplify the solving process and reduce the computational efforts. However, the DA model has limitations in estimating light propagation close to the sources and boundary or in the medium that contains the low-scattering region. Thus, more and more interest is now gradually turning to the light transport models based on solving the complete FD-RTE. Ren et al.[20] have employed the FD-RTE in the OT and have proved that the reconstruction algorithm based on the FD-RTE can reduce the cross-talk between the absorption and scattering coefficients significantly. Kim and Charette[21] have reported a sensitivity function-based conjugate gradient method to reconstruct the optical parameters with the FD-RTE, and proved that the frequency-domain data are more informative and better than steady data in separating internal parameters. Tarvainen et al.[26] have employed the Gauss-Newton reconstruction method in OT and obtained the reconstruction results of the case, including low-scattering regions based on the FD-RTE for the first time. Balima et al. have employed the discontinuous Galerkin finite element method[27] in OT to solve the FD-RTE, which can adapt the reconstruction algorithm to the complex geometry medium. Yamamoto and Sakamoto[28] have used a Monte Carlo perturbation calculation method to solve the FD-RTE in OT.

Although the complete FD-RTE has been used successfully in reconstructing optical parameters, the calculation of the gradient of the objective function with respect to the optical parameters is still problematic because almost all reconstruction algorithms are gradient-based methods in OT. The accuracies and efficiencies of these algorithms depend significantly on their gradient calculation techniques. In the past few years efforts have been made to calculate the gradient of the objective function, and many variant gradient-based methods have been developed to improve their performances in the sense of convergence speed and searching capability.[21,27,28] For instance, Lagrangian formalism, which can obtain the gradient via calculating the Lagrangian multipliers, has been widely used.[21,27] Additionally, the perturbation Monte Carlo technique[28] and gradient tree[29] have also been used to reconstruct the optical parameters. Although these algorithms used in OT can obtain the reconstructed images of absorption and scattering coefficients, the reconstructed results cannot show accurately the locations and shapes of the inclusions in the medium. This is due primarily to the fact that the obtained gradient lacks sufficient accuracy. To overcome this disadvantage, a more accurate reconstruction algorithm needs to be developed. The adjoint differentiation (AD) technique, which is always used in time-domain OT,[30,31] can calculate the gradient of the objective function accurately and efficiently without truncation error. However, the application of conventional AD technique requires the use of the numerical structure based on the time-step scheme, which is not contained in the frequency-domain model.

The motivation of the present study is to develop an accurate and efficient reconstruction algorithm based on AD technique to reconstruct the distributions of the absorption and scattering coefficients of inhomogeneous medium by using the FD-RTE. The remainder of this research is organized as follows. First, a finite volume method (FVM) formulation of the FD-RTE is well established and used as the forward model. Then, the objective function with the regularization term of the generalized Gaussian Markov random field (GGMRF) model is given to deal with the ill-posed problem. The multi-start conjugate gradient (MCG) method, which can speed up the convergence and enhance the accuracy of the reconstructed results, is employed to solve the inverse problem. A modified AD scheme which correlates the measured intensity with the inner optical parameters by using the collimated radiation is employed to calculate the gradient of the objective function. Finally, the simulation results of the forward model and the reconstruction results of the absorption and scattering coefficients in the inhomogeneous medium are discussed.

2. Forward model
2.1. Radiative transfer in frequency-domain

The FD-RTE is used as a forward model to simulate the propagation of photons in the participating medium. The FD-RTE in the direction Ω takes the following form:[32]

where , r denotes the location and ω is the modulation frequency c is the light speed in the medium, I (r,Ω,ω) is the radiative intensity at the position r in the direction Ω, μa (r), and μs (r) are the absorption and scattering coefficients, respectively, Φ (Ω′,Ω) is the scattering phase function that adopts the Henyey–Greenstein form in the present research:

with g being the anisotropic factor and θ being the scattering angle between the incident direction Ω′ and scattering direction Ω.

The intensity I within the medium is composed of the collimated intensity Ic and the diffuse intensity Id. The collimated component Ic obeys the following extinction law:

Substituting Eq. (3) into Eq. (1) yields

where Sc(r, Ω, ω) is the source term resulting from the collimated radiation and is given by

The collimated radiation Ic (r, Ωc, ω) can be expressed as

where tp is the pulse width and x is the distance of the collimated radiation in the direction Ωc.

The measurable quantity used for OT is the exitance on the boundary, thus

where n represents the normal vector of the boundary.

2.2. Solution method

The first and most essential stage in solving the complicated inverse radiation problem is to develop an accurate and efficient forward solution method. The FVM is employed to solve Eq. (4) in the present research for its suitability to any type of grid structure. The basic idea of FVM is to divide the computational domain into a large number of control volumes and ensure the radiative energy balances within each solid angle of the control volume. Integrating Eq. (4) in the control volume and control solid angle ΔΩm yields

where

is the source term and β = μa + μs is the extinction coefficient.

The left side of Eq. (8) can be written as

Converting the volume integral into a surface integral, the first term on the right side of Eq. (9) can be written as

where the Ac is the surface of the control volume VP, nj is the normal vector of Ac and Ωm represents the direction of the solid angle Ωm.

Assuming that the radiative intensities on the Ac are uniform in each Ωm and are equal to the intensity in the center of Ac, then

where represents the radiative intensity on the surface Ac,j in the solid angle Ωm. is the direction weight.

Representing the radiative intensity within the control volume VP by the intensity at the center point P of the VP in each solid angle Ωm, equation (8) can be expressed as

In each angle Ωm, the cell-surface intensity can be related to the cell-center intensity as

Then equation (13) can be written in a matrix form:[33]

where the J = E,W,S,N,F,B represents the center point of the control volumes, which neighbors the control volume VP. , , and can be written as

3. Inverse model

Theoretically speaking, the reconstruction of optical parameters is an inverse radiative problem that involves the estimation of the optical parameters inside the medium by measuring the exitance intensity on the boundary of the medium. The progress of reconstruction essentially makes the estimated exitance the measured exitance by changing the estimated optical parameters gradually. The inverse problem can be solved by finding the minimum value of the objective function.

3.1. Objective function

The objective function can be defined as a normalized-squared error between the predicted and measured exitance as follows:

where μ denotes the retrieval optical parameter and Ps,d (μ) and Ms,d are the estimated and measured values of the d-th detector with the s-th incident light source, respectively. To overcome the ill-posed nature of the inverse radiative problem, the GGMRF model[34] is adopted as the regularization term in the present research:

where p is the image sharpness parameter; σ is the scale parameter; N is the set of all neighboring pixel pairs; r is the neighboring position of s; bs−r is the weight coefficient. In this work, the CGMRF with p = 1.1 is used as the model which can provide good edge-preservation.[34] is used for the nearest neighboring position and is selected for the diagonal neighboring position. The value of σ is set according to the maximum likelihood (ML) estimate.[35]

3.2. MCG method

The reconstruction of optical parameters is equivalent to the minimization of the objective function. As a simple and effective optimization technique, the CG method is used to solve the minimization problem by iteration. During the iterations, the updating directional vectors of the optical parameters are determined by the gradient of the objective function, which can be written as

where dk is the descent direction vector and is the length of step. This procedure is repeated until the value of ∥f (μk+1 − F(μk)∥ is smaller than a predefined value ϕ.

The convergence speed of the CG method usually slows down with an increase of iteration times k. In particular, when μk is near the optimal solution, the convergence speed becomes intolerable. To address this issue, a multi-start iteration technique has been developed by our group, in which the descent direction vector dk becomes non-conjugated after the n-th iteration and d = −∇F (μ) have been used as the descent direction periodically to restart the iteration of the CG method.[36] This multi-start iteration technique, which stops the iteration at a fixed iteration number K and utilizes the presently reconstructed results as an initial guess to renew the iteration, is employed to improve the convergence efficiency of CG in solving the time-domain OT problem.[30] A flowchart of the MCG method is shown in Fig. 1. This method is employed to solve the frequency-domain OT problem in the present research.

Fig. 1. Flowchart of the MCG method.[30]
3.3. Calculation of the gradient

For the application of the MCG method, the core problem is to obtain the gradient dF(μ)/dμ of the objective function F(μ) with respect to the optical parameters μ. Given that the objective function is an implicit function of the absorption and scattering coefficients, the analytical expressions of the gradient dF(μ)/dμ are not available. Using the chain rule, the gradient of the objective function with respect to the optical parameters can be expressed as

where Ns and Ni represent the quantities of light sources and radiative intensities, respectively. The last term ∂F/∂ μ on the right side is zero, because the objective function F is not an explicit function of the optical property μ. The derivative ∂Ii/∂μ is the sensitivity of Ii to μ and can be obtained by partially differentiating the forward model with respect to μ. Equation (14) can be written as

Differentiating Eq. (14) with respect to μ yields

where the intensity Ii can be obtained by solving the forward model. ∂A/∂μ and ∂B/∂μ can be calculated easily using Eqs. (15a)–(15c).

For the time-domain model, the derivative can be calculated via the AD technique using the chain rule based on the time step scheme, and given by Eq. (23):[30]

where n is the time step.

Given that the objective function F is only the explicit function of the radiative intensity at the measured positions, the relation between the F and the inner radiative intensity must be established. The core significance of Eq. (23) is that it brings the information about inner intensity into the calculation of derivative dF/dIi by the derivative . In this way, the gradient of F(μ) with respect to the inner parameter can be obtained by Eq. (20). However, this time step scheme cannot be used for the frequency-domain model, and Eq. (23) is simplified into

where wm represents the weight associated with discrete direction vector Ωm and ξm is the cosine of the angle between n and Ωm.

Obviously, the derivative dF/dIi calculated by Eq. (24) is only a derivative with respect to the intensity at the measured positions. Thus, the gradient dF(μ)/dμ with respect to the inner parameter can be obtained only if the correlation between the measured intensity and the inner parameter is established. Notice that the matrix B in Eq. (22) contains the source term Sc if the collimated intensity Ic is included in the measured exitance. By Eqs. (5) and (6), the scheme of derivative dSc/du can be obtained as follows:

The derivative dSc/du contains the information about the inner optical parameters in the direction Ωc. Thus, the correlation between the measured intensity and inner parameter can be established by using Eq. (22) if each measured exitance contains the collimated intensity Ic. The gradient dF(μ)/dμ with respect to the inner parameter can then be obtained.

4. Results and discussion
4.1. Verification of the forward model

A one-dimensional uniform and isotropic slab medium (see Fig. 2) with a normally incident collimated irradiation on the left wall is considered to test the forward model. The thickness of the slab is L = 1.0 m, the extinction coefficient is β = 1.0 m−1, and the albedo is ω = 0.5. The non-dimensional pulse width of the incident radiation is . The slab medium is divided into 100 control volumes with 40 solid angles.

Fig. 2. Schematic diagram of the frequency-domain radiative transfer in a slab.

As shown in Fig. 3, the simulation results using FVM are compared with the results given by Liu and Hsu,[37] through using the discontinuous finite element method (DFEM). The simulation results of both methods are in excellent agreement with each other, which proves that the forward model employed in the present work can predict accurate results.

Fig. 3. The comparison between FVM and DFEM for (a) the hemispherical reflectance and (b) the hemispherical transmittance.

Next a two-dimensional uniform medium (see Fig. 4) with a collimated irradiation at the center of the left wall is considered to test the grid independence of the forward model. The size of the medium is 0.02 m × 0.02 m, the scattering coefficient is μs = 0.5 m−1, and the absorption coefficient is μa = 0.005 m−1. The pulse width of the incident radiation is tp = 1.67 ns.

Fig. 4. The schematic diagram of the frequency-domain radiative transfer in two-dimensional medium.

First, the number of solid angles is kept constant, and the FD-RTE is simulated by the forward model with grids of 11 × 11, 21 × 21, 31 × 31, and 41 × 41. Next, the number of grids is kept constant, and the FD-RTE is simulated with different quantites of the azimuthal and polar angles, respectively. The results (see Fig. 5) indicate that the transmittances with different numbers of grids and angles are in significant agreement with each other. Considering the fact that the forward model needs high computational efficiency to reduce the consumption of computation, the forward model with 80 solid angles and 21 × 21 grids is adopted in all the reconstruction tests.

Fig. 5. The verification of grid-independence for (a) the hemispherical transmittances with different grids. (b) The hemispherical transmittances with different azimuthal angles and (c) the hemispherical transmittances with different polar angles.
4.2. Description of the test model

The performance of the reconstruction algorithm is tested for retrieving the different distributions of the absorption and scattering coefficients of the participating medium with the anisotropic factor of g = 0.9. The frequency-domain measured exitance of the infrared light passing through a heterogeneous medium with inclusions are simulated with the forward FVM method. The detector readings are generated on the boundary of the medium. As shown in Fig. 6, the vertical incident point sources are set to be on the left and top wall of the phantom model in turn, and 19 × 4 detectors are distributed evenly over the boundary. The dense setting of the incident point sources aims to make sure that the measured exitance of each detector contains the incident intensity Ic. According to Ref. [38] the best quality of the reconstructed images can be achieved when a source-modulation frequency range of 400 MHz–800 MHz is used. Thus, the modulation frequency of the collimated source is set to be 600 MHz in the present research. The forward calculation model deals with an area of 2 cm × 2 cm with 80 discrete ordinate quadratures and 21 × 21 grids. The constant coefficients of the background medium are adopted as an initial guess. To test the stability and reliability of the reconstruction algorithm, the absorption and scattering coefficients of the media with different inclusions and noise data are reconstructed. The average time for one iteration step of a reconstruction on a PC with a 3.4 GHz i7-4770 CPU is about 20 min.

Fig. 6. The schematic diagram of the model with fibers.

To compare the overall accuracy of the reconstructed images the normalized root-mean square error (NRMSE) is introduced and defined as[39]

where is the retrieval value of the optical parameter at the r-th pixel, is the true value of the optical parameter at the r-th pixel, and N is the total number of pixels in the medium.

4.3. Reconstruction with different inclusions

To test the reliability of the proposed algorithm, the inclusions in the heterogeneous medium are set to be different locations and shapes. The absorption and scattering coefficients of the background medium are assumed to be μa = 5.0 m−1 and μs = 10.0 m−1. The absorption and scattering coefficients in the high extinction region are μa = 10.0 m−1 and μs = 15.0 m−1, and the absorption and scattering coefficients in the low extinction region are μa = 1.0 m−1 and μs = 5.0 m−1. The media with two square inclusions, two oval inclusions and ring inclusions are investigated in Test 1, Test 2, and Test 3 respectively.

The convergence process of the objective function in Test 1 is given in Fig. 7. The objective function decreases rapidly in the beginning of the iteration and reaches to a small value within 30 steps. So the reconstruction algorithm has high convergence efficiency for reconstructing the optical parameters based on the FD-RTE.

Fig. 7. Convergence process of the objective function.

As shown in Figs. 8–10, the distribution of the inclusions can be displayed accurately and clearly in the reconstructed images of the absorption and scattering coefficients. Both the locations and shapes of the different inclusions in the inhomogeneous medium are all reconstructed very well. As shown in Table 1, the values of NRMSE in the three tests are very small, which indicates that the proposed algorithm can reconstruct accurately the distributions of the absorption and scattering coefficients. The NRMSEs of the scattering coefficient are smaller than those of the absorption coefficient in all the tests, which indicates that the reconstructed images of the former have smaller errors than those of the latter. This can be attributed in part to the fact that the absorption coefficient is less than the scattering coefficient, and the influence of the absorption effect on the measured signals is less than that of the scattering effect.

Fig. 8. Real distributions of (a) absorption and (c) scattering coefficients and the reconstructed distributions of (b) absorption and (d) scattering coefficients for Test 1.
Fig. 9. Real distributions of (a) absorption and (c) scattering coefficients and reconstructed distributions of (b) absorption and (d) scattering coefficients for Test 2.
Fig. 10. Real distributions of (a) absorption and (c) scattering coefficients and reconstructed distributions of (b) absorption and (d) scattering coefficients for Test 3.
Table 1.

NRMSEs of the reconstructed images in three tests.

.
4.4. Reconstruction with measured errors

Measurement error is inevitable in the actual experimental measurement. To illustrate the performance of the proposed reconstruction algorithm, the random error added in the measured data is considered. The measured value with random error is obtained by adding the normal distributed error to the exact measured value:

where Ms,d and are the measured value and the exact measured value of the d-th detector, respectively, when the sth light source is working. The randn is a standard normal distributed random number. The standard deviation of the measured value σs,d, for a measurement error of γ% at 99% confidence, is determined as

where 2.576 arises from the fact that 99% of a normally distributed population is contained within ±2.576 standard deviation of the mean.

To test the robustness of the algorithm, 3%, 5%, and 10% random errors are added in the measurement data of Test 1, separately. The values of NRMSE (see Table 2) increase with an increase in the measurement error, which indicates that the qualities of the reconstructed images worsen when the measurement error increases. The NRMSE with measurement error keeps a relatively small value, and the reconstructed images (see Fig. 11) can display the shapes and locations clearly, thereby proving the good robustness of the algorithm proposed in this study.

Table 2.

The NRMSEs of the reconstructed images with considering measurement errors.

.
Fig. 11. Reconstruction results of absorption coefficient (top row) and scattering coefficient (bottom row) with measurements errors of 3% errors (left column), 5% (middle column) and 10% (right column).
5. Conclusions

An accurate reconstruction algorithm for OT based on the FD-RTE is developed in the present research. The forward model using the FVM formulation of the FD-RTE is employed to predict detector readings on the surface of a medium, with the given absorption and scattering coefficients. An objective function with a regularization term is defined by the predicted detector readings and simulated measurement detector readings. The GGMRF model is used as the regularization scheme to overcome the ill-posed nature of the inverse problem. The MCG method is employed as an optimization technique to reconstruct the absorption and scattering coefficients. A modified AD scheme using the collimated intensity is employed to calculate the gradient of the objective function with respect to optical parameters. Considering the test results, the following conclusions can be drawn.

Reference
1Chen X LLei Y Z 2015 Chin. Phys. 24 030301
2Cui MGao XZhang J 2012 Int. J. Therm. Sci. 58 113
3Sun S CQi HZhao F ZRuan L MLi B X 2016 Appl. Therm. Eng. 98 1104
4Wu X CJiang H YWu Y CSong JGrehan GSaengkaew SChen L HGao XCen K F 2014 Opt. Lett. 39 638
5Zhang BQi HRen Y TSun S CRuan L M 2014 J. Quantum Spectrosc. Radiat. Transfer 133 351
6Wu F YYang X B 2012 Chin. Phys. 21 057803
7Zirak A RKhademi M 2007 Opt. Commun. 279 273
8Liu Z MLiu W QGao M GTong J JZhang T SXu LWei X L 2008 Chin. Phys. 17 4184
9Lohmann A W 1991 Opt. Commun. 86 365
10Musculus M P BPickett L M 2005 Combust. Flame. 141 371
11Wyatt J SCope MDelpy D TRichardson C EEdwards A DWray SReynolds E O1990J. Appl. Phys.681086
12Klose A DHielscher A H 2002 J. Quantum Spectrosc. Radiat. Transfer 72 715
13Wang F QTan J YShuai YTan H PChu S X 2014 Int. J. Heat Mass Tranfer 78 7
14Gassan S AAndreas H H 2003 J. Electron. Imaging 12 594
15Elaloufi RCarminati RGreffet J J 2002 J. Opt. A.: Pure Appl. Opt. 4 S103
16Quan HGuo Z 2004 Opt. Express 12 449
17Boulanger JCharette A 2005 J. Quantum Spectrosc. Radiat. Transfer 91 297
18Zirak A RKhademi M 2007 Opt. Commun. 279 273
19Qi HQiao Y BSun S CYao Y CRuan L M 2015 Math. Probl. Eng. 2015 1
20Ren KBal GHielscher A H 2006 SIAMJ Sci. Comput. 28 1463
21Kim H KCharette A 2007 J. Quantum Spectrosc. Radiat. Transfer 104 24
22Joshi ARasmussen J CSevick-Muraca E MWareing T AMcGhee J 2008 Phys. Med. Biol. 58 2069
23Chu MVishwanath KKlose A DDehghani H 2009 Phys. Med. Biol. 54 2493
24Pogue BTestorf MMcBride TOsterberg UPaulsen K 1997 Opt. Express 1 391
25Gibson A PHebden J CArridge S R 2005 Phys. Med. Biol. 50 R1
26Tarvainen TVauhkonen MArridge S R 2008 J. Quantum Spectrosc. Radiat. Transfer 109 2767
27Balima OFavennec YBoulanger JCharette A 2012 J. Quantum Spectrosc. Radiat. Transfer 113 805
28Yamamoto TSakamoto H 2016 Opt. Commun. 364 165
29Meng JWang J JHuang X WLiu R J 2006 J. Quantum Spectrosc. Radiat. Transfer 102 181
30Qiao YQi BChen H QRuan L MTan H P 2015 Opt. Commun. 351 75
31Klose A DHielscher A H 1999 Med. Phys. 26 1698
32Arridge S R 1999 Inverse Probl. 15 R41
33Mishra S CChugh PKumar PMitra K 2006 Int. J. Heat Mass Transfer 49 1820
34Saquib S SHanson K MCunningham G S 1997 Proc. SPIE 3034 369
35Saquib S SBouman CSaue K A 1998 IEEE Trans. Image Proc. 7 1029
36Yuan XSun W1997Theory and Methods of OptimizationBeijingScience Press191
37Liu L HHsu P F 2008 J. Quantum Spectrosc. Radiat. Transfer 105 357
38Kim H KNetz U JBeuthan JHielscher A H2009SPIE BiOS: Biomedical Optics. International Society for Optics and Photonics 71742B-710.1117/12.809409
39Ye J CWebb K JBouman C AMillane R P 1999 J. Opt. Soc. Am. 16 2400